\documentclass{ctexart}
\usepackage{standalone}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{tikz}
\usepackage{listings}
\usepackage{amssymb}


\usetikzlibrary{decorations.markings}

\title{step-3}


\author{金郁香 \\信息与计算科学 3220103054}

\begin{document}

\maketitle

\section{有限元方法的基本设定}
这是第一个我们真正使用有限元来计算的例子。我们将解决一个简单的泊松方程,其边界值为零
\begin{align}
  -\Delta u=f\quad\quad\quad\quad \mathrm{in~}\Omega,
  \\u=0\quad\quad\quad\quad \mathrm{on~}\partial\Omega.
\end{align}
我们将在正方形区域$\Omega=[−1,1]^2$上求解此方程，如果你已经了解了有限元方法的基础知识，你将记得我们需要采取以下步骤来用有限维近似来逼近解u。具体来说，我们首先需要导出上述方程的弱格式，方法是将方程左乘以一个测试函数$\varphi$，并在域$\Omega$上进行积分（我们将在下面回来解释为什么要从左边乘而不是从右边乘）:
\begin{align}-\int_{\Omega}\varphi\Delta u=\int_{\Omega}\varphi f.\end{align}

这可以通过分部积分得到:
\begin{align}\int_\Omega\nabla\varphi\cdot\nabla u-\int_{\partial\Omega}\varphi\mathbf{n}\cdot\nabla u=\int_\Omega\varphi f.\end{align}

测试函数$\varphi$必须满足与问题边界条件相同的条件（数学上来说: 它需要来自我们寻求解的集合的切空间），因此在边界上$\varphi = 0$，因此我们要求的弱格式为
\begin{align}
(\nabla\varphi,\nabla u)=(\varphi,f),\end{align}

其中我们使用了常见的符号表示 $(a,b)=\int_{\Omega}a\mathrm{~}b$.然后，该问题要求一个函数u，使得对于所有适当空间中的测试函数$\varphi$，上述语句都成立（这里是空间$H^{1}$）。

当然，在一般情况下，我们无法在计算机上找到这样的函数，而是寻求一个近似解$u_h(\mathbf{x})=\sum_jU_j\varphi_j(\mathbf{x})$,其中$U_j$是我们需要确定的未知扩展系数（该问题的“自由度”），而$\varphi_i(x)$是我们将使用的有限元形状函数。为了定义这些形状函数，我们需要以下内容：

1.网格：需要定义一个网格，用于确定形状函数。你已经在 step-1 和 step-2 中学习了如何生成和操作描述网格的对象。

2.有限元：需要选择一个有限元来描述形状函数，该有限元针对参考单元（在deal.II中通常是单位区间[0,1]、单位正方形$[0,1]^{2}$或单位立方体$[0,1]^3$）定义形状函数。在 step-2 中，我们已经使用了类型为 FE\_Q<2> 的对象，它表示常规的拉格朗日元素，通过插值在支撑点上定义形状函数。其中最简单的一个是 FE\_Q<2>(1)，它使用一阶多项式。在二维中，它们通常被称为双线性函数，因为它们在参考单元的两个坐标中都是线性的（在一维中，它们是线性的，在三维中，它们是三线性的）。（然而，在 deal.II 的文档中，我们经常不区分这一点，而是统称这些函数为“线性函数”。）

3.DoFHandler：需要使用一个 DoFHandler 对象，在网格上枚举所有自由度，以参考单元描述有限元对象提供的基础。你已经在 step-2 中学习了如何做到这一点。

4.映射：需要告诉如何从参考单元上有限元类定义的形状函数得到实际单元上的形状函数。默认情况下，除非你明确指定，deal.II将使用一个（双线性、三线性）映射来完成这一步骤，因此在大多数情况下，你不需要担心这一步。

通过这些步骤，我们现在有了一组函数$\varphi_i$，并且我们可以定义离散问题的弱形式：找到一个函数$u_h$，即找到上面提到的扩展系数$U_j$，使得
\begin{align}
(\nabla\varphi_i,\nabla u_h)=(\varphi_i,f),\quad\quad\quad i=0\ldots N-1.\end{align}

请注意，在这里我们遵循从零开始计数的约定，这在C和C++中很常见。这个方程可以重写为一个线性系统，如果你插入表示$u_h(\mathbf{x})=\sum_jU_j\varphi_j(\mathbf{x})$，然后观察到

\begin{align}
(\nabla\varphi_i,\nabla u_h)& =\left(\nabla\varphi_i,\nabla\left[\sum_jU_j\varphi_j\right]\right)  \\
&=\sum_j\left(\nabla\varphi_i,\nabla\left[U_j\varphi_j\right]\right) \\
&=\sum_j\left(\nabla\varphi_i,\nabla\varphi_j\right)U_j.
\end{align}

因此，问题可以表述为：找到一个向量U，使得
\begin{align}AU=F,\end{align}

其中矩阵A和右手边F定义为
\begin{align}A_{ij}=(\nabla\varphi_i,\nabla\varphi_j),\\F_i=(\varphi_i,f).\end{align}


\section{组装矩阵和右手边的向量}

现在我们知道我们需要什么（即：保存矩阵和向量的对象，以及计算$A_{ij}和F_i$的方法），我们可以看看如何实现：

我们需要一种方式来计算积分。在有限元方法中，最常用的方法是使用数值积分（quadrature），即将积分替换为对每个单元上一组数值积分点的加权求和。首先，我们将整个域$\Omega$上的积分分解为每个单元上的积分：

\begin{align}A_{ij}=(\nabla\varphi_i,\nabla\varphi_j)=\sum_{K\in\mathbb{T}}\int_K\nabla\varphi_i\cdot\nabla\varphi_j,\\F_i=(\varphi_i,f)=\sum_{K\in\mathbb{T}}\int_K\varphi_if,\end{align}

然后通过数值积分来近似每个单元的贡献：

\begin{align}
&A_{ij}^K =\int_K\nabla\varphi_i\cdot\nabla\varphi_j\approx\sum_q\nabla\varphi_i(\mathbf{x}_q^K)\cdot\nabla\varphi_j(\mathbf{x}_q^K)w_q^K,  \\
&F_{i}^{K} =\int_K\varphi_if\approx\sum_q\varphi_i(\mathbf{x}_q^K)f(\mathbf{x}_q^K)w_q^K, 
\end{align}

其中$\mathbb{T}\approx\Omega$ 是表示域的一个三角剖分，$\mathbf{x}_q^K$是第q个积分点在单元K上的位置，$\mathbf{w}_q^K$是第q个积分权重。进行这个过程需要不同的部分，我们将逐个讨论它们。

首先，我们需要一种方式来描述数值积分点的位置$\mathbf{x}_q^K$和权重$\mathbf{w}_q^K$。通常情况下，它们通过与形状函数相同的方式从参考单元映射，即使用MappingQ1类隐式描述，或者如果您明确指定，通过Mapping派生类中的一个类进行映射。参考单元上的位置和权重由继承自Quadrature基类的对象进行描述。通常，人们选择数值积分公式（即一组点和权重），使得数值积分与矩阵中的积分完全相等；这可以通过高斯积分公式实现，该公式在QGauss类中实现。

然后，我们需要一些方法来帮助我们在单元K上评估$\varphi_i(\mathbf{x}_q^K)$。这就是FEValues类所做的：它接受一个有限元对象来描述参考单元上的$\varphi$，一个数值积分对象来描述数值积分点和权重，以及一个映射对象（或隐式地使用MappingQ1类），并在位于K上的数值积分点上提供形状函数及其导数的实际单元K上的值，以及用于积分的其他所需信息。

\section{求解线性系统}
我们在这里要做的是采用1952年的一个想法：共轭梯度法（Conjugate Gradient method），或简称为 "CG"。CG是一个"迭代"求解器，因为它形成了一个向量序列，可以收敛到精确解；事实上，在没有舍入误差的情况下，经过N次这样的迭代，如果矩阵是对称和正定的，它就可以找到精确解。

这些迭代求解器的一个重要组成部分是,我们需要指定要解决线性系统的容差-实质上,是关于我们对愿意接受近似解的误差的一个声明。线性系统Ax=b的精确解x的近似解$\tilde{x}$的误差定义为$\|x-\tilde{x}\|$,但这是一个我们无法计算的量,因为我们不知道精确解x。相反,我们通常把残差,定义为$\|b-A\tilde{x}\|=\|A(x-\tilde{x})\|$,作为一个计算措施。然后我们让迭代求解器计算越来越精确的解$\tilde{x}$,直到$\|A(x-\tilde{x})\|\leq\tau$。一个实际的问题是$\tau$应该取什么值。在大多数应用中,设为
\begin{align}&\tau=10^{-6}\|b\|\end{align}
  是一个合理的选择。
  
\section{关于实现}
下面这个程序显示了大多数有限元程序的基本结构:
\begin{lstlisting}[language=C++, frame=single]
class Step3
{
  public:
    Step3();
    void run();

  private:
    void make\_grid();
    void setup\_system();
    void assemble\_system();
    void solve();
    void output\_results() const;

    Triangulation<2>     triangulation;
    FE\_Q<2>              fe;
    DoFHandler<2>        dof\_handler;

    SparsityPattern      sparsity\_pattern;
    SparseMatrix<double> system\_matrix;
    Vector<double>       solution;
    Vector<double>       system\_rhs;
};
\end{lstlisting}


\section{关于类型}
deal.II通过命名空间types中的别名定义了一些整数(integral)类型。(在上一句话中，"积分(integral)"一词被用作形容词，与名词integer相对应。它不应该与代表曲线或曲面下的面积或体积的名词"integral" 相混淆。形容词integral在C++世界中被广泛使用，如 "积分类型"(integral type)、"积分常数"(integral constant)等。）特别是，在这个程序中，你会在几个地方看到types::global\_dof\_index：一个整数类型，用来表示自由度的全局索引，即在定义在triangulation之上的DoFHandler对象中特定自由度的索引（相对于特定单元中特定自由度的索引）。对于当前的程序（以及几乎所有的教程程序），你将有几千个到几百万个全局未知数（对于元素，在2d的每个单元上有4个未知数，在3d有8个未知数）。因此，允许为全局DoF指数存储足够大的数字的数据类型是无符号整数unsigned int，因为它允许存储0到略高于40亿的数字（在大多数系统中，整数是32位的）。事实上，这就是types::global\_dof\_index的内容。
  
 using types::global\_dof\_index = typedef unsigned int
  
\section{描述计算过程和具体参数}
Step-3是一个求解二维拉普拉斯方程的程序。它使用deal.II库进行有限元数值计算。程序中定义了一个名为`Step3`的类，其中包含了求解步骤的各个函数，包括创建网格、设置系统、装配系统、求解和输出结果。

问题的描述：

1. 创建网格：使用'GridGenerator::hyper\_cube()'函数生成一个单位超立方体网格，并通过'triangulation.refine\_global()'函数进行全局细化。

2. 设置系统：使用'DoFHandler'类分配自由度，并通过'DynamicSparsityPattern'类创建稀疏矩阵模式。

3. 装配系统：使用高斯积分公式对每个单元进行数值积分，计算单元刚度矩阵和右手边向量，并将其累加到整个系统矩阵和右手边向量中。

4. 求解：使用共轭梯度法（CG方法）对线性方程组进行求解，得到数值解。

5. 输出结果：使用DataOut类将数值解写入文件。

具体参数：

- 建立的网格是一个二维的单位超立方体网格。

- 有限元使用的是二次Lagrange多项式（FE\_Q<2>）。

- 使用高斯积分公式进行数值积分，积分精度与有限元阶数相同。

- 边界条件使用零函数。

- 求解器为共轭梯度法（CG方法），迭代次数上限为1000，收敛精度为1e-12。

- 输出的结果保存在名为"solution.eps"的文件中。 

\section{生成图片}

\includegraphics[width=1\textwidth]{solution.pdf}

\end{document}
